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Abstract. High resolution reconstruction of complicated objects from incom¬ 
plete and noisy data can be achieved by solving modulation equations iteratively 
under physical constraints. This direct demodulation method is a powerful tech¬ 
nique for dealing with inverse problem in general case. Spectral and image restora¬ 
tions and computerized tomography are only particular cases of general demod¬ 
ulation. It is possible to reconstruct an object in higher dimensional space from 
observations by a simple lower dimensional instrument through direct demodula¬ 
tion. Our simulations show that wide field and high resolution images of space 
hard X-rays and soft 7 -rays can be obtained by a collimated non-position-sensitive 
detector without coded aperture masks. 


1. Introduction 


An observation to an object in space x with an intensity distribution f{x) can be seen as 
a modulation process of signal from object by observation instrument, the observed data 
d{u) is an output after the modulation, where to denotes the parameters determining the 
state of observation. The modulation equation relating the data to object can be written 
as 


p{(jo,x)f{x)dx = d{u) 


( 1 ) 


where a value p{uj, x) of the integral kernel (modulation function) is the modulation 
coefficient or response coefficient of the instrument to a point x of the object space 
during an observation cj. 


Generally, the data space cj is not the same kind of the object space x. Eor 
an imaging formation system consisting of position-sensitive detectors, the data space 
cj = a;' is a two-dimensional space correspondent to the object space x and the observed 
data d{x') is an image of the object f{x), if the modulation function (point-spread 
function) is space-invariant, p{x', x) = p{x' — x), then the modulation equation (image 
formation equation) can be written as the following convolution formula 


J p{x' — x)f{x)dx =p*f = d{x') . 


( 2 ) 


The image formation Equation (2) is a particular case of the general modulation Equation 
(I), and the image restoration (restoring the object f{x) from the observed image d{x')), 
then, is a particular case of object reconstruction. 

Dividing object space into N bins, for M observed values d{k), k = 1,...,M, 
the modulation Equation (1) or image formation Equation (2) in discretization form 
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Figure 1. 


Projection along line L in two-dimensional Radon transform. 


constitute an algebraic equation system 

N 

'^p{k,i)f{i) = d{k) {k = l,...,M). (3) 

i=l 

In this paper we introduce the direct demodulation method of reconstructing object 
f from observed data d. This method can be applied to any kind of observation described 
by Equation (3). For a start we will give a brief summation of various reconstruction 
methods in classihcation. 


2. Methods of Reconstruction 
2.1. Linear Transformation 

For observations by an image formation system, taking Fourier transforms of both sides 
of the image formation Equations (2) and taking advantages of the convolution theory 
of the Fourier transform, we have J-p ■ Tf = Ed, then with the inverse transform the 
object can be derived as / = E~^{EdlEp). 

In the case of observed values d being projections, i.e. line integrals of object f{x) 
along lines L : a = ^ • x (see Figure 1) 

f f{x)dx= j f{x)S{a — ^-x)dx = d{a,(l)) 

then d is a sample of the Radon transform (Radon 1917), d = Hf, with the inverse trans¬ 
form we can reconstruct the object from projections, / = lZ~^d. The Radon transform 
is the mathematical basis of the computerized tomography (CT) technique. 

The cross-correlation function 

c{x) = J p{u,x)d{u)du (4) 

is a distribution in the object space. In general case we can use the correlation calcula¬ 
tion to perform transformation from data space to object space. Correspondent to the 
modulation Equation (3), the discrete correlation transform is 

c{i) = ^p{k,i)d{k) . 

k 


( 5 ) 
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There exists another form of correlation transform defined as 

A ^ 

-P)id{k)-d) 

k=i 

= A'[Mj2p(.k,i)d{k)-J2piku)J2d{k)]- ( 6 ) 

k k k 

If the object is an isolated point source with intensity fs at point ig, then the correlation 
function c*{i) = fs ■ Cs{v.,is)-, where Cs{v,is) is the resolution function for point is 

Cs{i]is) = A'[M^p{k,i)p{kAs)-^p{k,i)^p{kAs)\ ■ (7) 

k k k 

The distribution of c*{i) will appear a bulge around the bin is, provided that the ob¬ 
servations get necessary information of object (the observations are properly arranged 
and the number M of observations is large enough). The normalization constant can 
be chosen as A' = I/[M'f2k(pik,i)Y ~ (^kP{k,i)y‘\, then the maximum value of the 
correlation function, c*{is), equals the source strength, fg. The width of the resolution 
function Cs{i) represents the resolution ability of the observations, which can be defined 
as the intrinsic resolution of the observations. 

2.2. Statistical Reconstruction 

Another important category of reconstruction methods is to use a certain statistic (e.g. 
the maximum likelihood ratio (Pollock et al. 1981)) instead of the object itself to rep¬ 
resent main characters of object by some degree or to use values of a statistic (e.g. the 
quantity) instead of the observation modulation equations in selecting possible object 
distributions. If the object distribution can be described by a theoretical or empiri¬ 
cal expression f{i) = g{i;a), where a = (ai,...,a/) are I undetermined parameters, the 
least-squares reconstruction g{i; a) {i = 1,..., N) satisfies the condition 

X^(a) = a)]/o-fc}^ = min 

k i 

with ak being the statistical error of d{k). The maximum-entropy method (for example, 
see Gull &: Daniell 1978) is to choose a solution /(i) {i = 1, ...,N) from all satisfying the 
statistical criterion 


M 

X^if) = '^{['^Pik,i)fii) - dik)]/ak}‘^ = M 

k=l i 


by the condition that the information entropy 

Si f) = - Ylfi fi = max. 

i 


2.3. Direct Demodulation 

The most straightforward way to evaluate an object / from observed data d should be 
directly solving the modulation equation connecting d with /, i.e. direct demodulation, 
which we will discuss in this paper. 

Using a statistic instead of original modulation equations to make reconstruction 
will cause loss of information containing in the equations about object and observations, 
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then that of sensitivity and resolution ability of reconstruction. The least-squares method 
is model dependent, not applicable in general case. Reconstruction by Fourier or Radon 
transform is mathematically based on modulation equation, but only can work on a 
special kind of observation, i.e. image formation or projection. In application of Fourier or 
Radon transform to incomplete and noisy data in practice, in particular to that with poor 
statistics and low ratio of signal to noise, in order to derive a reasonable reconstruction 
a certain treatment (e.g. frequency space filtering) should be applied, causing loss of 
information. Correlation calculation can perform transformation from data space to 
object space in general case, thus being a useful tool in dealing with reconstruction 
problems; but the correlation function c{i) or c*{i) is not a good representation of an 
object, there may exist regions with negative c*{i) even if the object is defined only in a 
positive space. 

Before introducing the direct method of reconstruction, we will firstly compare 
various methods of solving a linear equation system. 


3. Solutions of Linear Algebraic Equations 

The modulation equation system (3) can be rewritten in matrix form as 

Pf = d, ( 8 ) 

where modulation matrix P = {p(/c,i)}, object vector / = data vector d = 

{d{k)}, i = 1,...,A, /c = 1,..., M. 

3.1. Mathematical Solution 

If the modulation matrix P is non-singular and M = N, the exact solution of Eq.(8), 
f = P~^d, can be found by Gaussian elimination or triangular resolution. Iterative 
method can also be used to solve Eq.(8). An iteration algorithm in common use is 
the Gauss-Seidel method with successive relaxation. The approximate solution for l-th. 
iteration can be calculated by one of the following two formulas 

i-l N 

a{di - + (1 - + (1 - a)/f ( 9 ) 

jr = l j=*+l 

i-l N 

— {di-^Pijff- ^ + {I- ( 10 ) 

J=1 3=i+l 

where the relaxation factor 0 < a < 1. It can be proved that the Gauss-Seidel iteration 
converges to the exact solution if the matrix P is symmetry and positive definite. 

The mathematical solution of a modulation equation system satisfy the equations 
accurately. But there always exist errors in modulation equations (measurement errors, 
statistical fluctuation and noises in data and errors in determining modulation coeffi¬ 
cients). A modulation equation system in a practical problem is usually a large set of 
linear equations, its mathematical solutions will seriously deviate from the true object 
distribution due to errors in the equations. A one-dimensional example is displayed in 
Figure 2; (a) is an object distribution /, a line at zq on a uniform background; (b) the 
ideal data dQ{k) = where the triangular distribution presents the shape of 

the modulation function p{k, zq ); ( c ) the simulated observed data d, a Monte Garlo sam¬ 
ple from (b); (d) and (e) are the correlation distribution c(i) and c*(z) respectively; (f) 
presents the mathematical solution of the modulation equations Pf = d by Gauss-Seidel 


fP 

Ai) 
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Figure 2. Solutions of modulation equations, (a) Object spectrum, (b) The¬ 
oretical observed spectrum, (c) Simulated observation data, (d) Correlation 
distribution c(i). (e) Correlation distribution c*{i). (f) Mathematical solution 
of modulation equations, (g) Least-squares solution of modulation equations, 
(h) Solution of correlated modulation equations of L = 2. (i) Iterative solution 
with nonnegative constraint. All horizontal axes are in the same scale, the hor¬ 
izontal line in each plot indicates the zero level, the vertical axes are bounded 
from —200 to 600 for (a), (e), (g), (h), (i), that from —100 to 300 for (b), (c), 
(d) and from —10000 to 10000 for (f). 


iterations. From Figure 2f we can see that the mathematical solution of the modulation 
equations with statistical errors is in an ill-conditioning feature with violent oscillation, 
including a lot of nonphysical values of / < 0. 

3.2. Statistical Solution 

If the number of observed values is greater than that of the undetermined parame¬ 
ters, M > N, the system (3) is over deter mined, we can use the least-squares criterion 
J2k['I2iP{ku)f{i) ~ d{k)]‘^ = min to derive the normal equations as 

N 

= ci{i') {i'= ( 11 ) 

i=l 

where pi{i', i) = J2kPik^ i')p{k, i), the right side of the normal equation is the correlation 
function ci(i') = c(i') = '^kP{k,i')d{k). The mathematical solution of the normal Equa¬ 
tions (11), i.e. the least-squares solution of the modulation Equations (3), can be taken 
as an estimation of object in a sense of least-squares (in a statistical sense). 

The normal equations are just the correlation transform of the original modulation 
Equations (3), we can call them correlated modulation equations. Eor the case oi M < N 



6 


TI-PEI LI AND MEI WU 


we can also derive equations of the same form by means of correlation transformation. 
The correlation transformation can be repeatedly performed (Li & Wu 1992, 1993). 
Starting from the original modulation eqs. (3), after L-fold correlation transforms we 
can obtain the L-order correlated modulation equations as 

N 

= CL{i') (/= l,...,iV), (12) 

i=l 

where PL(i',i) = Y.jPL-i{jD')VL-i{jD)NL{i') = EjPL-i(j)*')cL-i(j)- 

Figure 2g and 2h show the solutions of correlated modulation equations of L = 1 
and L = 2, respectively (statistical solutions of modulation equations). 


3.3. Physical Solution (Constrained Iterative Solntion) 


Setting reasonable physical constraints in the iterative process of solving a modulation 
equation system (3) or correlated modulation equation system (12) will help us obtain 
an estimation of the object distribution satisfying some necessary physical conditions (Li 
& Wu 1992,1993). Corresponding to a certain object, we can set up a lower bounds b 
and upper bounds u, for any approximate solution 

if < b{i), let f^^\i) = b{i), 

if f^’'\i) > u{i), let f^^\i) = u{i). (13) 


The final solution should be normalized by the condition J2k'l2iP{l^D)f{i) = Z)fc^(^)- 
For many objects, nonnegativity is a necessary constraint, i.e. b{i) = 0, (i = 1,...,A^). 
Figure 2i is a constrained iterative solution for data (c) under the nonnegative condition 
after 10^ iterates. The constrained iteration process has a good property of convergence 
and is not sensitive to small errors in the modulation function. We found the approximate 
solution after 10^ or more iterates or that starting from completely different initial values 
had not any recognized difference with Figure 2i. When the modulation function used 
in iteration being widened in a direction and narrowed in the other direction by 20% 
separately, no obvious change was found in the constrained iterative solution, and when 
the whole modulation function being widened by 20% also almost the same hgure of 
solution was obtained but just the peak height increased a little. 

Exerting the constraints of lower and upper bounds is a kind of nonlinear control 
in the iterative precess, proper nonlinear control can depress the influence of noise, fluc¬ 
tuation and other errors in data effectively and strengthen convergence of the iterative 
process, helping to derive a satisfactory reconstruction in a physical sense. There exist 
various iterative algorithms applicable to seeking physical solutions. The modulation 
matrix of the correlated modulation Equations (12) is symmetry and positive defi¬ 
nite, satisfying convergence condition of the successive relaxation Gauss-Seidel iteration. 
Another useful algorithm is the Richardson-Lucy iteration (Richardson 1972; Lucy 1974) 
based on the Bayes theory 


/«(,) = 

k 


p{k, i)d{k) 




(14) 


Solutions from probability iterations satisfy the nonnegative condition and nor¬ 
malization condition of total flux automatically. However, only nonnegative constraint is 
not enough to control the iterative process to produce a satisfactory reconstruction from 
many observations, in particular from that of poor statistics and low ratio of signal to 
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noise for complicated objects. In the case of object space including negative area, the 
simple nonnegative constraint can not be used at all. 

Relaxation functions to modify the correction term gradually as iteration proceeds 
were suggested for spectral deconvolution (for example, see Jansson 1984). In oder to 
fulfil a need for setting various kind of physical constraints necessary in practice we 
proposed to set boundary constraints directly on each approximate solution as shown by 
expression (13) and found that a direct demodulation method in the general case can be 
constituted by this technique, which we will introduce in the following section. 


4. Direct Demodulation Method 

It is obvious that in solving modulation Equations (3) or correlated modulation Equations 
(12) taking the intensities of the diffuse background in the object space as lower 
bounds, i.e. 

h = fb (15) 

can make more precise control than using simple nonnegative constraint, giving a better 
estimation of the object. 

In the case of no priori knowledge of the diffuse background, it is needed to estimate 
it from observed data. Eirst we use a successive procedure similar to CLEAN (Hogbom 
1974) to subtract all contribution of apparent discrete sources with intensities greater 
than fm from the observed data d and derive background data db-, where fm can be 
taken as the minimum intensity of discrete sources to be reconstructed. Solving the 
modulation equations for db, Pfb = db, by an iterative method we finally get fb- In 

iterations we can simply smooth each approximate solution or use some kind of 
continual condition like 

uii) = f^\i - 1) + <5, b(t) = fPii -l)-5 (16) 

to depress sudden change of intensity at any point i from its close neighbour one i — 1 
in the object space, where 5 = fm/nh with rib being the linear size of minimum diffuse 
structure reconstructed by number of bins. 

Subtracting contribution of apparent discrete sources from data can be performed 
with the aid of cross-correlation technique. Eind out the maximum point ig of the correla¬ 
tion map c*{i) {i = 1,..., N) for data d and the intensity excess fg above the background 
level at ig, remove the contribution of fg from the data by db{k) = d{k) — p{k,ig)fg, 
{k = 1,...,M). The procedure is repeated for db till fg < fm- The excess fg can be 
determined by the condition — [c*(*0 “ fsCsii',is)])‘^ = min, where the sum 

for i' is accumulated over a region around ig within the width of the resolution function, 
Cb{i') is the extrapolated value of correlated intensities c*{i) outside the region to the 
point i!. 

An example of restorations of a line source superimposed on an inclined background 
is shown in Figure 3; (a) is the object distribution, (b) the simulated data, (c) Richardson- 
Lusy iterative restoration (i.e. a constrained iterative solution of modulation equations 
with only the nonnegative condition), (d) restoration by the direct demodulation for 
correlated modulation equations of L = 1. 

For many reconstruction problems in practice, performing iterations in two steps 
with the constraint condition (16) and (15) separately by the Gauss-Seidel formula with 
successive relaxation to the correlated modulation Equations (12) of L = 1, or by the 
Richardson-Lucy formula to the original modulation Equations (3), one can get satis¬ 
factory reconstructions. In the case of data having good statistics, solving the original 
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Figure 3. Direct demodulation, (a) Object spectrum, (b) Simulated data, 
(c) Richardson-Lusy restoration, (d) Direct restoration. 




Collimator 


Scintillator 


PMT 


Figure 4. A collimated scintillation telescope.Direct demodulation. 

modulation Equations (3) by constrained Gauss-Seidel iterations will give high resolution 
reconstruction; but in the other hand, for observations with severe noise and fluctuation 
and bad intrinsic resolution, higher order correlated modulation Equations (12) can serve 
to get better reconstruction. In iterations for reconstructing diffuse background corre¬ 
lated modulation equations of higher order are also recommended. 


5. Modulation Imaging 

The direct demodulation method can be applied to dealing with various kind of re¬ 
construction problems, i.e. to solving the inverse problems of the general modulation 
Equations (3), where the data space determined by apparatus and states of observa¬ 
tion can be completely different with the observed object space x, and the number of 
observations, M, is not necessary to be equal to or greater than that of undetermined 
parameters, N (Li et al. 1993). A general demodulation technique make it possible to 
design an experiment by simple instrument to realize the same objective as more com¬ 
plicated one does, for example, to obtain a two-dimensional image by a one-dimensional 
non-position-sensitive detector ( non-imaging system). 
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Figure 5. Modulation imaging, (a) True intensity distribution, (b) Cross¬ 
correlation map from simulated data obtained by a collimated telescope shown 
in Fig.4. (c) Reconstruction by direct demodulation, (d) Direct reconstruction 
for the case of variable background. 


We made Monte Carlo simulations to study the imaging capability of a simple colli¬ 
mated hard X-ray telescope consisting of two 700 cm^ non-position-sensitive scintillators, 
each in viewed by just a single photomultiplier tube (PMT) and modulated by a slat col¬ 
limator of 5° X 0.5'^ and 0.5" x 5° (FWHM) aperture separately (sketched in Figure 4). 
The object scene (Figure 5a) in our simulations is the galactic center region containing 
eight point sources which was observed by a coded mask telescope flown on the Spacelab 
2 mission (Skinner et al. 1987). The strongest source intensity was 7.05 x 10“^ cm“^ 
s“^, the ratio of the strongest source intensity to the weakest was 100:1, the background 
was taken as 0.089 counts cm“^ s“^ as estimated from the observation. A 24 hours 
survey by the collimated telescope was simulated: pointing to various sky points over a 
6° X 6 ° region around the Galactic Center separated by 1° along the galactic longitude 
I or latitude b between two successive ones, during each pointing producing the output 
counts of each detector for eight different places of the telescope xy plane rotated around 
the z axis by Monte Carlo simulation. Figure 5b is the cross-correlation map from the 
simulated data. A reconstruction by the direct demodulation technique is shown in Fig¬ 
ure 5c, which seems even better than the simulated image of the same object scene by 
a 1473 cm^ coded mask telescope with 93 x 99 detector pixels for 24 hours observation 
(Lei et al. 1991). 

For studying the effect of background variation on modulation imaging, we doubled 
the background intensity when the telescope pointing to {l,b) = (1.0°,—1.0°) but no 
recognizable change was produced in the resultant image. In addition we assumed that 
starting from the pointing {l,b) = (1.0°,—1.0°) the background intensity was linearly 
increasing with the passage of time up to three times and then gradually decreasing 
back to the initial level, the whole variation lasted one hour, the corresponding direct 
reconstruction is shown in Figure 5d. 

The modulation imaging has a capability to simultaneously reconstruct both ex¬ 
tended and discrete features in object. We simulated scan observations of the object 
scene of 10° x 10° shown in Figure 6a by the collimated telescope (Figure 4) with step 
of 1°, Figure 6b is the direct demodulation image from the simulated data. 
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Figure 6. (a) Object scene, (b) Direct reconstruction from simulated data 

obtained by a collimated telescope shown in Figure 4. 


6. Discussion 


The direct demodulation, directly solving modulation equations and directly controlling 
each iterative output by physical conditions, is a flexible technique; in accordance with 
different characters of object, instrument and state of observation and different require¬ 
ment for reconstruction corresponding algorithms can be designed, and various kind of 
priori knowledges can be easily included in modulation equations or constraint conditions 
as well. 

The modulation Equation (1), observed data d{u) being an weighted sum (inte¬ 
gral) of an object f{x) over any region with a weight function (integral kernel) x) of 
any form, can represent a measurement process in very general case, thus the direct de¬ 
modulation method can be applied to general reconstruction problems. The CT imaging 
is only a special case of general modulation imaging; the observed data are projections 
of object; in other words, the CT technique can be used only to a special case of the 
modulation Equation (1) where the integral kernel (modulation function) is a kind of 6 
function. 

Each iterative calculation in direct demodulation is directly based on the modu¬ 
lation equation. Obviously the direct approach can use the information containing in 
the modulation equations more sufficiently than any indirect one through a statistical 
criterion or cross-correlation transformation, that is the main reason of the resolution of a 
direct reconstruction being able to be much better than the intrinsic resolution (compare 
Figure 2i with 2e and Figure 5c with 5b). For ideal data without error and with neces¬ 
sary information of object, i.e. the modulation matrix P being non-singular (observations 
properly arranged) and its rank M > N (number of observations large enough), solving 
the modulation equations will precisely reconstruct the object, e.g. solving Pf = 
can exactly restore Figure 2a from Figure 2b no mater how wide the resolution function 
p{k,io) is. 

The main difficult of the direct method is the illness of solution causing by errors in 
data. Introducing physical constraints to make nonlinear control in the process of solving 
modulation equations and using correlation transformation properly can depress the effect 
of errors effectively and then make high resolution direct reconstruction realizable. 

The physical solution resulted from solving a modulation equation system under 
physical constraints is neither a solution in a strictly mathematical sense nor the statis¬ 
tical solution in a least-squares sense. Neither the modulation equations nor the normal 
equations are satisfied by the physical solution. Studying the mathematical properties 
of physical solutions is a difficult problem of nonlinear mathematics. The quality of a 
direct reconstruction (sensitivity, resolution, uncertainty, etc.) depends on the quality 
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and quantity of information about the object containing in the observed data as well as 
the ratio of signal to noise, and relates to the ability of the signal modulation pattern de¬ 
pressing noises. The generality of the direct demodulation technique implies that various 
kind of apparatus (detectors and modulators) constituted by different principles can be 
selected for observing the same object. Until now no general mathematical method can 
guide us to designing modulators and arranging observations, even not a clear standard 
to judge the quality of a reconstruction. How many observations are necessary for recon¬ 
struction depends on properties of observation system, arrangement of observations and 
characteristics of observed object (size of reconstructed space and complexity of object), 
and also relates to requirements for reconstruction. Just from resolution consideration 
we may arrange observations properly and take the number M of observations suffi¬ 
ciently large to make the intrinsic resolution function (7) fine enough (comparing with 
the achievable limit of the instrument). 

Designing detector and modulator and estimating statistical errors and other un¬ 
certainties in a reconstruction can be aided by simulation calculations. Statistical uncer¬ 
tainties can be evaluated by the bootstrap method (Efron 1979; Diaconis & Efron 1983), 
estimating other uncertainties and comparing different observation instruments and ar¬ 
rangements can be made by Monte Carlo simulations. In performing direct demodulation 
by a computer, the computation is quite fast as the algorithm is constituted by only sim¬ 
ple arithmetic operations and conditional control statements, although a larger space 
of the main storage is often needed to store modulation coefficients and reconstruction 
results. It is feasible in many practical cases to study the direct demodulation through 
simulation calculations for Monte Carlo samples or bootstrap samples. 

The authors thank Prof. Jun Nishimura for helpful discussion. This work was 
supported by the National Natural Science Eoundation of China. 
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